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Abstract 

Goal, Scope, and Background. As Life Cycle Assessment (LCA) 
and Input-Output Analysis (IOA) systems increase in size, com¬ 
putation times and memory usage can increase rapidly. The use 
of efficient methods of solution allows the use of a wide range 
of analysis techniques. Some techniques, such as Monte-Carlo 
Analysis, may be limited if computational times are too slow. 
Discussion of Methods. In this article, I describe algorithms that 
substantially reduce computation times and memory usage for 
solving LCA and IOA systems and performing Monte-Carlo 
analysis. The algorithms are based on well-established iterative 
methods of solving linear systems and exploit the power series 
expansion of the Leontief inverse. The algorithms are further 
enhanced by using sparse matrix algebra. 

Results and Discussion. The algorithms presented in this article 
reduce computational time and memory usage by orders of mag¬ 
nitude, while still retaining a high degree of accuracy. Lor a 
3225x3225 LCA system, the algorithm reduced computation 
time from 70s to 0.06s while retaining an accuracy of 10- 3 %. 
Storage was reduced from 166 megabytes to 1.8 megabytes. The 
algorithm was used to perform a Monte-Carlo analysis on the 
same system with 1,000 samples in 90s. I also discuss various 
issues of power series convergence for general LCA and IOA 
systems and show that convergence will generally hold due to 
the mathematical structure of LCA and IOA systems. 
Conclusions. By exploiting the mathematical structure of LCA 
and IOA iterative techniques substantially reduced the compu¬ 
tational times required for solving LCA and IOA systems and 
for performing Monte-Carlo simulations. This allows more wide¬ 
spread implementation analysis techniques, such as Monte-Carlo 
analysis, in LCA and IOA. 

Recommendations and Perspectives. It is suggested that algo¬ 
rithms, such as the ones described in this article, should be im¬ 
plemented in LCA packages. Various checks can be used to verify 
that computational errors are kept to a minimum. 

Keywords: Algorithms; input-output analysis (IOA); matrix in¬ 
verse; Monte-Carlo; numerical approaches; power series expan¬ 
sion; sensitivity analysis 


Introduction 

Life Cycle Assessment (LCA) and Input-Output Analysis 
(IOA) are becoming two prominent methodologies for 
analyzing the environmental impacts of consumer activities 
and production processes. LCA has traditionally focussed 


on detailed studies of particular products and services (Bau¬ 
mann and Tillman, 2004); such as paper cups (Hocking 
1991). In contrast, IOA has traditionally focussed on more 
aggregated products and economic activities; such as house¬ 
hold consumption (Hertwich 2005). Unfortunately LCA 
often suffers from ill-defined system boundaries (Lave et al. 
1995) and consequently IOA is often used to increase the 
systems boundaries of LCA studies (Suh et al. 2004). In¬ 
creasingly, LCA studies are either supplemented with IO data 
(e.g. Suh 2004) or performed almost entirely using IO data 
(Hendrickson et al. 1998, Joshi 1999). 

With the increased use of LCA and IOA, there is also an 
increase in the size of the underlying databases. Currently, 
the biggest IO databases contained about 500 industry sec¬ 
tors (Bureau of Economic Analysis 2005) and some com¬ 
mercial LCA databases have data on around 2,500 proc¬ 
esses (Frischknecht and Jungblush 2004). Given this, it is 
not uncommon for hybrid LCA systems to require the solu¬ 
tion of equations with 2,000 to 3,000 variables. Even with 
modern computers, computational times for systems of this 
size can limit the applicability of some analysis techniques. 
For instance, detailed studies of parameter uncertainty us¬ 
ing Monte-Carlo analysis may require the repeated solution 
of a system thousands of times (Huijbregts et al. 2003). For 
the wide spread use of numerical approaches, such as Monte- 
Carlo analysis, it is essential that computational resources 
are kept to a minimum. 

Currently, there are two broad techniques for solving LCA 
systems (Heijungs and Suh 2002); the sequential approach 
and the matrix approach. The sequential approach starts 
with one central process and linearly scales the required prod¬ 
ucts and processes; this is repeated further upstream. The 
matrix approach converts the LCA into a system of linear 
equations which are solved using standard linear algebra 
(Heijungs 1994, Heijungs and Suh 2002). IOA systems are 
only solved using matrix approaches (Miller and Blair 1985). 
Consequently, the use of the matrix approach allows the 
similarities of LCA and IOA to be exploited in hybrid LCA 
studies (Joshi 1999, Nakamura and Kondo 2002, Suh 2004, 
Peters and Hertwich 2006, Suh 2006). Further, by express¬ 
ing LCA in the same matrix form as IOA allows easy appli¬ 
cation of the many established tools developed for IOA (e.g. 
Miller and Blair 1985, Lenzen 2003). For these reasons the 
remainder of this article will focus only on the matrix ap¬ 
proach for LCA. 

Generally, solving an LCA system using the matrix approach 
requires the solution of a system of linear equations (Heijungs 
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1994, Heijungs and Suh 2002). Direct methods of solving a 
system of linear equations include matrix inversion and 
Gaussian elimination (Nicholson 1990). However, it is well 
established that these direct methods are computationally 
inefficient (c.f. Watkins 2002). Slow computation times have 
made some LCA interpretation methods less viable. For in¬ 
stance, Heijungs et al. (2005) limit the number of samples 
used for Monte-Carlo analysis due to computational con¬ 
straints as opposed to accuracy considerations. It is well es¬ 
tablished that for large systems of linear equations iterative 
methods are considerably faster then direct methods of so¬ 
lution (Watkins 2002). 

This article discusses algorithms for the efficient solution of 
LCA and IOA systems and demonstrates the methods using 
Monte-Carlo analysis. Instead of using direct methods to solve 
systems of linear equations, I apply well established iterative 
methods which are computationally superior without com¬ 
promising accuracy. Further, I discuss how the structure of 
IOA and LCA equations ensure the convergence of iterative 
methods. Most of the analysis draws upon well established 
results in the numerical analysis of linear systems (e.g. Berger 
and Saibel 1957, Stewart 1973, Seneta 1973, Hackbusch 
1994, Datta 1995, Saad 1996, Varga 2000, Watkins 2002). 
I demonstrate the algorithms using Monte-Carlo analysis. 

1 Solution of Linear Systems 

For a given LCA, a system of linear equations can be writ¬ 
ten in the same form as in IOA 1 

x = Ax + y (1) 

where each element of the vector x is the output of the vari¬ 
ous production activities, each element of the vector y is the 
demand of products on the LCA (the functional unit), and A 
is the normalized technology matrix. In the normalized form 
the z-th column of A represents the inputs of products required 
to produce one unit of the z-th product. I assume that there is 
the same number of production activities denoted by n. 1 The 
normalized form in (1) is subtly different to the un-normal- 
ized equations used by some authors (e.g. Heijungs and Suh 
2002). As I will show, to use iterative methods of solution as 
described below it is essential to have a normalized system for 
convergence to hold in an arbitrary system 3 . In Appendix A 
I describe how to construct a normalized system of linear 
equations from an un-normalized system. 

Analytically the solution of (1) can be obtained by matrix 
inversion, 

x = (I - A) -1 y = Ly (2) 


1 A detailed description on how to construct an LCA system as a system of 
linear equations can be found in Heijungs and Suh_(2002); however, 
they construct an un-normalized system in the form As = y. This is dis¬ 
cussed further in Appendix A. 

2 For LCA, Heijungs and Suh (2002) describe how to deal with systems 
with different numbers of processes and products. In IOA the make and 
use system describes how to deal with secondary production (Miller and 
Blair 1985, United Nations 1993). 

3 Heijungs and Suh (2002), pp. 108-109, give an illustrative example of this. 


In the Leontief inverse, L = (I-A)- 1 , the column i represents 
the required products to produce a unit of product z. The 
Leontief inverse has received relatively little attention in LCA, 
despite its importance in IOA (Miller and Blair 1985). If the 
environmental impacts per unit output, F, are known for 
each activity, then the total environmental impacts for the 
functional unit, y, are 

f=F(I-A)~iy (3) 

where the rows of F represent different environmental 
stressors and the columns represent the different produc¬ 
tion processes. Consequently, f is a column vector with each 
element representing a different environmental stressor. 

While calculating the inverse, as in (2), or using Gaussian 
elimination are effective direct means of solving (1) they are 
not computationally efficient (Watkins 2002). Rather, itera¬ 
tive methods are often used to solve (1). Interestingly, ex¬ 
pressing (2) in the form (1) is the starting point for many 
iterative numerical methods 4 of solving equations of the form 
(2) (c.f. Berger and Saibel 1957, Stewart 1973, Seneta 1973, 
Hackbusch 1994, Datta 1995, Saad 1996, Varga 2000, Wat¬ 
kins 2002). These methods are generally of the form 


QO 

X = Yj X i < 5 ) 

z=o 

Solving this iterative equation for z —> °° leads to the solution 

x = y + Ay + A 2 y + A 3 y + ... (6) 

This power series expansion is commonly used in IOA as a 
way of interpreting the Leontief inverse (Miller and Blair 
1985, United Nations 1999). The z-th term in the expansion 
represents the required products in the z-th production layer. 
In the zero-th production layer, the output y is from the 
manufacturing site or sites. To produce y an input of Ay is 
required into the first tier production layer. To produce Ay 
an input of A(Ay) = A 2 y is required into the second tier pro¬ 
duction layer, and so on through the infinite expansion. 

In terms of IOA, the power series expansion 

(I-A)- l = I + A + A 2 + A 3 *... (7) 

was also derived independently by Waugh (1950) and it was 
further shown that the series converges if A is a non-nega¬ 
tive matrix and if the row sum of A is less than one, 

n 

Zl41 <1 fora11 / w 

i=l 


4 Many different iterative methods are defined based on how the iterative 
'A' matrix is constructed from the starting point x = Ly. For simplicity I 
will assume that a system in the form (2) is rewritten as (1). The refer¬ 
ences describe the many methods of constructing an 'A' matrix. 
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This convergence condition holds for IOA systems in mon¬ 
etary units, but is not directly useful in mixed unit IOA or 
LCA since the matrices generally do not satisfy these condi¬ 
tions. A more useful result from linear algebra is that (7), 
hence (6), converge when the spectral radius 

p(A) = {max (|d|) for all eigenvalues, A, of A} (9) 

is less than one 5 (e.g. Seneta 1973, Datta 1995, Varga 2000). 
Further, this result does not require A to be non-negative. 
The spectral radius can also be used as an alternative to the 
Hawkins-Simon conditions (Hawkins and Simon 1949) for 
testing the existence and uniqueness of positive solutions to 
(1) (Wood and O'Neill 2002). 

The iterative scheme (4) and power series expansion (6) offer 
an alternative way to solve (1). For most LCA systems it is 
generally straightforward to verify convergence (see below), 
however, of more interest is whether (4) is a computationally 
superior way of solving (1), particularly for larger systems. 

Computational methods for solving linear systems of equa¬ 
tions have been studied in considerable detail in the past. 
For an arbitrary system of equations in the form (1) it is 
generally found that direct calculation of the inverse takes 
| n 3 + 0(n 2 ) computational flops to solve 6 , Gaussian elimi¬ 
nation takes f n 3 + 0(n 2 ) flops, and iterative methods usu¬ 
ally take Cn flops for some constant C which depends on 
the number of iterations required for convergence 7 (Hack- 
busch 1994, Saad 1996, Watkins 2002). The constant C can 
be large and so iterative methods are generally better for 
large n. If the linear system of equations is sparse 8 then these 
computational times are found to be significantly reduced 
(Hackbusch 1994, Saad 1996, Watkins 2002), although the 
computational times depend on the degree of sparsity. Over¬ 
all, it is generally found that on most large problems the 
savings in computation times and storage achieved by using 
iterative methods 'are truly spectacular' (Watkins 2002). 

2 Algorithms 

In this section I describe simple algorithms that apply the 
iterative scheme in (4) and then compare the computational 
times to direct methods of solution. I also discuss the issues 
of convergence and accuracy. 

Using a simple recursive procedure it is possible to imple¬ 
ment (4) as in Algorithm 1. To reduce computational time 


5 In words, the spectral radius is the largest absolute value of the 
eigenvalues of A. In physical terms, if all the eigenvalues of A lie inside 
a disc of radius one then the power series converges. This is analogous 
to A- converging for Ixl < 1. A number A is called an eigenvalue of A if 
it satisfies Av= Avtor some non-zero vector v. Eigenvalues and eigen¬ 
vectors have numerous applications in linear algebra (c.f. Nicholson 1990). 

6 The convention is usually to assume that the computational work is the 
same for additions, subtractions, multiplications, and divisions. The big- 
O notation, Ofn 1 ), essentially states that, in the limit of large n, the com¬ 
putational time of the algorithm varies proportional to n 1 . 

1 It is also possible to determine the required iterations in advance (Seneta 
1973, Datta 1995, Varga 2000); although a more efficient method is to 
add terms until an acceptable error or number of iterations is reached 
(see later). 

8 A sparse matrix has a high proportion of zero elements compared to 

non-zero elements. 


and memory usage, the properties of sparse matrices are 
exploited 9 (Hackbusch 1994, Saad 1996). Algorithm 1 takes 
inputs of A, y, the desired accuracy of the solution ( toler¬ 
ance ), and the maximum number of iterations to perform 
(max-iterations). The first four lines of the algorithm ini¬ 
tialize the variables; test is given a value to ensure that at 
least one iteration occurs. The algorithm then recursively 
constructs the output until the desired accuracy is reached 
or the maximum iterations are exceeded. The index vari¬ 
able i is the current tier of the expansion (6), x, holds the 
value of the current tier which is constructed recursively from 
the previous tier, see (4) and (6), x contains the total output, 
and the accuracy is measured relative to the previous itera¬ 
tion (test). It is also possible to calculate the Leontief inverse 
by replacing y with the identify matrix in Algorithm 1. 

Typically the objective of an LCA is to calculate the envi¬ 
ronmental impact for a given functional unit, y. This allows 
minor modification of Algorithm 1 to allow direct calcula¬ 
tion of environmental impacts, Algorithm 2. The structure 
of the algorithm is the same as Algorithm 1. The main dif¬ 
ference is that an extra step is used to calculate the environ¬ 
mental impact. In terms of computational time, Algorithm 2 
is superior due to the faster computation of test , and fur¬ 
ther, it can be used to give both the output and the environ¬ 
mental impacts. 

Algorithm 1: Calculate_Olitpijt(A, y, tolerance , maxjterations) 
i<- 0 
x i <— y 

X <— Xj 

test <— tolerance 

while (test > tolerance) and (i < maxjterations) 
i <— i+1 
x i A *x n 

X <— X + X/ 


IMI 

end 

return x 

Algorithm 2: CALOJLATE_lMPAcr(F, A, y, tolerance , maxjterations) 

l <-0 

Xi<~y 

f <—F * x 

test <— tolerance 

while (test > tolerance) and (i < maxjterations) 
i <— i+1 
x, ^ A * x,_, 

X <— X + X; 

f t <—F * x 

test 1 

|/| 

end 

return f t 

9 Many computational packages, for instance Matlab, have optimized rou¬ 
tines for sparse matrix calculations. 
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2.1 Computational times 

To give an indication of the computational times required 
for direct methods of solution compared to indirect itera¬ 
tive methods of solution, several calculations were performed 
using Matlab Version 7 operating under Microsoft Windows 
XP (Table 1). Of course, depending on the machine, the 
operating system, the operating state of the machine, and 
the code used for the calculations, these times will vary. The 
A matrix for this calculation was 3225x3225 and 1.4% of 
the matrix elements were non-zero. The data set is a hybrid 
LCA system based on ecoinvent (Stromman et al. 2006). 
The functional unit, y, has one non-zero element. The table 
clearly illustrates that the power series expansion gives more 
than sufficient accuracy and computational times three or¬ 
ders of magnitude faster than conventional techniques. The 
error in Table 1 compares the direct method of solution with 
the iterative method of solution. For a matrix system of this 
size it is not feasible to analytically calculate the solution 
and so the direct method is also a numerical approximation 
(see below for more details). The table shows that the toler¬ 
ance gives a good indication of the error of the iterative 
method compared to the direct method. 

In terms of memory usage, the use of sparse matrix algebra 
saves considerable space. For the system used here with 64- 
bit (8 byte) real numbers requires 3225x3225x8 = 83 mega¬ 
bytes of storage. The same storage would be required for 
the identity matrix. If the matrices are stored as sparse ma¬ 
trices then A requires 1.8 megabytes and the identity matrix 
requires 51 kilobytes. These savings may also lead to re¬ 
duced computational time since system resources are avail¬ 
able for other computational activities. 

2.2 Convergence 

In general, the convergence of iterative schemes can be veri¬ 
fied using a variety of techniques; however the structure of 
LCA and IOA systems makes convergence easier to verify. 
With some conditions it is possible to prove that for an arbi¬ 
trary LCA or IOA system p(A) < 1 and convergence holds. 
For the case of IOA monetary units it has been shown that 
converge always holds (Waugh 1950). In general, LCA and 
IOA systems are in mixed units so the row and column sums 


cannot be used to ensure convergence as they may be greater 
than one. In the case of mixed units convergence follows if 
A is a positive matrix. Given the structure of (1), it follows 
that I and A are a regular splitting of the system (e.g. Saad 
1996, pp 107) and thus if L> 0 then it follows that p(A) < 1 
(Saad 1996, also see O'Neil and Wood 1999, Wood and 
O'Neil 2002). 

In some situations A may have non-negative coefficients. 
This usually arises due to assumptions on waste-streams and 
secondary production in both LCA (Heijungs and Suh 2002) 
and in IOA (Miller and Blair 1985). This generally occurs in 
very few elements and experience has shown this is gener¬ 
ally not a problem in IOA (United Nations 1999). From a 
theoretical point of view, it may be difficult to guarantee 
convergence if elements are arbitrarily non-negative. Al¬ 
though, from a practical point of view it would be expected 
that convergence holds. If convergence did not hold, then 
the interpretation of the power series expansion (6) would 
not hold which would suggest that the A matrix was not 
constructed correctly. That is, the divergence of A' as i —» °° 
implies that in the z-th production layer requires an infinite 
input of products. Thus it is likely that p(A) < 1 holds for all 
realistic LCA and IOA systems. From this practical point of 
view, if convergence did not hold then the validity of A would 
be questionable. 

Given possible doubts about the convergence of A it is possi¬ 
ble to use a variety of methods to test convergence (Saad 
1996). For instance, it could be numerically verified that 
p(A) < 1, the matrix norm is less than one, or the matrix is 
diagonally dominant. Again, direct calculation is computa¬ 
tionally slow and iterative methods are usually used, such 
as the power method 10 (the algorithms can be found in Datta 
1995, Wood and O'Neil 2004). Using the power method, 
together with the properties of sparse matrices, reduces com¬ 
putation times of the spectral radius by orders of magnitude 
whilst still retaining accuracy (Table 2). Thus, computa¬ 
tionally it is straightforward to verify if p(A) < l. n 


10 Not to be confused with the power series expansion. 

11 Although, the pragmatic may skip this step and if (6) does not converge 
then it probably means that r. 


Table 1: The computational times for CalculateJmpact. The asterisk, *, indicates the number of iterations to get an error of ICh 3 percent. The direct 
calculations were performed using the Matlab functions inv which performs a matrix inverse and \ which uses Gaussian elimination. The error is relative 
to the direct method and the tolerance is from Algorithm 2 and compares the last two iterations 


Tiers 

Time 

Time 

Error (%) 

Tolerance (%) 

full (s) 

sparse (s) 



Direct, inv 

70 

40 

- 

- 

Direct, \ 

14 

4.3 

- 

- 

10 

0.6 

0.03 

2.3 

1.4 

20 

1.2 

0.05 

3.3x1 O' 2 

1.8x10 -2 

27* 

1.5 

0.06 

2.0x10“ 3 

9.7x10" 4 

50 

2.9 

0.10 

2.7x10“ 7 

1.3x10 -7 

100 

5.0 

0.15 

3.7x10“ 9 

8.6x10“" 
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Table 2: The computational times for determining the spectral radius. The asterisk, *, indicates the number of iterations to get an error of 10- 3 percent. The 
Direct calculation was performed using the Matlab call max(abs(eig(A))). The error is compared to the direct solution and the tolerance is the comparison 
of the last two iterations. The power method was used for the iterative solutions (e.g. Datta 1995, Wood and O'Neil 2004) 


Iterations 

Time 

Time 

Error (%) 

Tolerance (%) 

non-sparse (s) 

sparse (s) 

Direct 

100 

- 

- 

- 

10 

0.6 

0.03 

40 

11 

20 

1.3 

0.06 

14 

11 

50 

3.2 

0.15 

8.4x10 -4 

4.1X10 -3 

51* 

3.2 

0.17 

3.7x10" 4 

4.7x10" 4 

100 

6.5 

0.31 

9.7x10 -8 

2.5x10“ 7 


2.3 Accuracy 

Another issue of importance is the accuracy of the method. 
I have designed the algorithms to continue until a desired 
accuracy or maximum number of iterations is reached. If 
the matrix is ill-conditioned then both the direct and itera¬ 
tive solution may be sensitive to small changes in the func¬ 
tional unit. The condition number of a matrix is given by 

^HWIWf < 10 > 

where |A| is a matrix norm and it follows that k(A)> 1. A 
matrix is defined as ill-conditioned if k(A) » 1 (Watkins 
2002). However, an ill-conditioned matrix is not always sen¬ 
sitive to small changes (Watkins 2002), particularly for sparse 
systems where the matrix columns consists mainly of zeros 
and are therefore nearly linearly-dependent. 

The condition number for the matrix used in this article is 
of the order 10 23 which suggests that the problem is highly 
ill-conditioned. However, by back-substituting x back into 
(1) gives an error of 10 _2 % showing that the result is in fact 
accurate. In addition, the results presented in Table 1 show 
that the system is converging to the same solution for differ¬ 
ent numbers of iterations. A number of tests were also per¬ 
formed by perturbing the functional unit and the system was 
found to be stable. 

The reason this system is stable despite being ill-conditioned 
is probably related to the fact that it is very sparse (1.4% of 
elements are non-zero) and consequently the columns are 
nearly linearly-dependent (causing the matrix to appear sin¬ 
gular). It was also found that the data directly from ecoinvent, 
which forms the foundation for the data used in this article, 
is ill-conditioned. Consequently, for these systems it is best 
to verify that the solution is accurate by back-substitution 
rather than relying on the condition-number; this approach 
is also computationally faster. 

3 Monte-Carlo Analysis 

Monte-Carlo analysis is a technique that propagates known 
parameter uncertainties through a calculation to give an 
uncertainty distribution on the output variables (Morgan 
and Henrion 1990). Consequently, Monte-Carlo analysis is 


an ideal method for quantifying parameter uncertainty in 
TCA studies (e.g. Huijbregts et al. 2001, 2003) and IOA 
studies (e.g. Bullard and Sebald 1988, Tenzen 2001). How¬ 
ever, it seems Monte-Carlo analysis is rarely performed in 
TCA and IOA studies; one possible reason for this is a belief 
that Monte-Carlo analysis is computationally slow (c.f. 
Heijungs et al. 2005). In this section, I apply Algorithm 3 to 
reduce the computational times of conventional Monte-Carlo 
studies. Ciroth et al. (2004) have a similar motivation, but 
do not use the matrix approach which makes there method 
not as viable for hybrid LCA systems. 

The algorithms presented in the previous section are par¬ 
ticularly relevant in performing a Monte-Carlo analysis. The 
accuracy of a Monte-Carlo simulation generally increases 
with the square root of the number of samples taken; a reli¬ 
able Monte-Carlo analysis might consist of 10,000 or more 
samples 12 . The algorithms presented above give the oppor¬ 
tunity for large increases in the sample size in a Monte-Carlo 
simulation. 

I applied a simple Monte-Carlo analysis as described in Al¬ 
gorithm 3. The algorithm has several functions which need 
explaining. The variable N is the number of samples in the 
Monte-Carlo calculation. I have assumed that F, A, and y 
are stored as sparse matrices. For performing loops on a 
sparse matrix it is considerably faster to loop through the 
array representation of a sparse matrix than the full matrix 
(e.g. Saad 1996). The function Find() is an operation that 
returns an array representation of a sparse matrix 13 (Saad 
1996); that is, the row and column indices and the value of 
the non-zero entries are returned. The function Sparse)) con¬ 
verts the array representation of a matrix back to the sparse 
matrix representation 14 . In the procedure Random)), the func¬ 
tion PDF( V, o) returns randomly distributed numbers with a 
given probability distribution; in the example here, a nor¬ 
mal distribution. The error is taken as three standard devia¬ 
tions which covers 99% of the distribution. 


12 The sample size will vary for different problems and it is difficult to make 
a generalization of how many samples are adequate. Heijungs et al. (2005) 
suggest 500 to 1,000 samples. This may not always be an adequate 
number of samples. Huijbregts et al. (2003) use 10,000 samples. 

13 As in the Matlab find routine. 

14 As in the Matlab sparse routine. 
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Algorithm 3: Monte_Carlo(P, A, y, N, error, tolerance, 
maxjterations) 

procedure Random(R, C, V, cr) 
comment: This procedure gives a probability 
distribution to the non-zero elements 
for e 1 to Length(R) 

V new [e] ^?DF(V[e], a[e]) 
end 

return Sparse (R, C, V new ) 

main 

throws’ A columns? ^values! FlND(A) 

rows’ Vcolumns’ Vvaluesl FlND( V) 

error/3 A values 
F a f- error/3 F values 

for / <— 1 to N 

A sample <- Random(A toim , A columns’ ^ values’ ^<7) 

^sample Random(F tows , F 

columns’ F values’ ^<j) 

E[i] <- CALCULATE_lMPACT(F Mmpfc , A sample , y, tolerance, 
maxjterations) 

end 

return E 

For the Monte-Carlo analysis I only apply distributions to 
the non-zero elements. The easiest way to operationalize this 
process is to convert the matrices into a sparse representa¬ 
tion. This is performed in the first two lines of Algorithm 3; 
for each non-zero element in A, A rows contains the row in¬ 
dex, A columns contains the column index, and A va[ues the value 
of the non-zero element. For the demonstrative purpose of 
this example, each element in the matrices is given a con¬ 
stant percentage error with a normal distribution 15 . The 
standard deviation of the coefficients is given by F a and A ^ 
the error is assumed to be three standard deviations. Once 
the data has been initialized, the algorithm then calculates 
the emissions for each sample of the Monte-Carlo. The pro¬ 
cedure Random produces a new sample distribution of the 
matrices and then Calculate_Impact is used to determine 
the environmental impacts with a given accuracy ( tolerance) 
for each sample. The algorithm returns a vector containing 
the environmental impacts for each Monte-Carlo sample. 

Algorithm 3 was demonstrated on the system described 
above and the results are shown in Table 3. I used a nor¬ 
mally distributed error of 25%, with the error defined as 


15 A log-normal (or other) distribution may be more appropriate in many 
practical studies(e.g. Huijbregts et al. 2003; Bullard and Sebald 1988). 


three standard deviations. The Monte-Carlo results stabi¬ 
lize at around 1,000 samples, suggesting that 1,000 samples 
are adequate for this system. This calculation takes approxi¬ 
mately 90 seconds for the 3225x3225 system described 
above. Although, the time penalty for using more samples is 
not great. With smaller systems it is likely that Monte-Carlo 
simulations for 10,000 or more samples will take only a 
matter of minutes. 

The computationally slowest part of the algorithm is in the 
calculation of the randomly distributed matrices, that is, the 
function Random. There are many options for speeding up 
this function. For instance, a smaller sample size can be used 
by dividing the sample space into equiprobable intervals; 
such as in Latin Hypercube Sampling 16 (Morgan and Henrion 
1990). Sampling techniques give a more uniform distribu¬ 
tion for a given sample size allowing the sample size to be 
reduced; however, they are also computationally expensive. 
Thus, sampling techniques are only used if the computation 
times with a reduced sample size offset the time to do the 
sampling. For our example problem, the procedure 
Calculate_Impact is sufficiently fast so that it is slower to 
reduce the sample size using sampling techniques. 

4 Conclusion 

In this article I have demonstrated the use of iterative algo¬ 
rithms for solving LCA and IOA systems. The algorithms 
were demonstrated on a sparse 3225x3225 hybrid LCA sys¬ 
tem and they performed several orders of magnitude faster 
than direct methods of solution. Several issues relating to 
the convergence of the algorithm were discussed and it was 
shown that the structure of LCA and IOA systems will gen¬ 
erally lead to convergence. The algorithm was then used to 
perform a Monte-Carlo simulation on the same system. It 
was demonstrated how efficient methods of solving the LCA 
system can be used to employ a large number of samples in 
a Monte-Carlo analysis. The algorithms presented here make 
it possible to perform Monte-Carlo analysis on most LCA 
and IOA systems in a matter of seconds to minutes. 


Acknowledgements. Anders Hammer Stromman provided the data 
for the simulations performed in the article. Two reviewers provided 
useful comments which greatly improved the article. 


16 The Matlab routine Ihsnorm uses Latin Hypercube Sampling. 


Table 3: The computational times and accuracies for the Monte-Carlo simulation. The mean values estimate Ffl-AJy, Std is the standard deviation, and 
the percentage error is given by three standard deviations 


Samples 

Mean 

Std 

Error (%) 

Time 

exact 

103.5 




10 

102.4 

7.1 

20.8 

00:00:01 

100 

104.2 

9.2 

26.4 

00:00:10 

1,000 

103.3 

8.4 

24.4 

00:01:30 

10,000 

103.5 

8.7 

25.1 

00:15:00 

100,000 

103.5 

8.7 

25.1 

02:30:00 
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Appendix A: Normalized systems of linear equations for LCA 

In this article I have used a normalized system of linear equa¬ 
tions to represent an LCA system. Some authors prefer to 
use an un-normalized system. In this appendix I briefly de¬ 
scribe the differences between the two systems. 

Consider an LCA system written in an un-normalized form 
(Heijungs and Suh 2002) 

As = y (11) 

where y is the demand on the LCA system, s is a scaling vector 
for the different processes, and A is the un-normalized 'tech¬ 
nology' of the LCA system. I have assumed that each process 
in the un-normalized system has only one output through al¬ 
location (Heijungs and Suh 2002) and the processes are or¬ 
dered to have the output on the diagonal of A. 

The normalized system, (1), can be obtained by using 

A = [d~A)D~' =I-AD~ l (12) 

where D is a matrix with the diagonal elements of A on the 
diagonal and hence the scaling vector s is replaced by the 
output, 

x = Ds (13) 

Further, the emission intensity also needs to be normalized, 
compare with (3), so that 

F = FD~ 1 ( 14 ) 

It is worth elaborating on a few differences in the normal¬ 
ized and un-normalized systems. In the un-normalized sys¬ 
tem, (11), the output of the processes are placed on the di¬ 
agonal of A and the off-diagonal terms are negative to 
represent 'inputs' into the production processes. The ele¬ 
ments of s represent the 'share' of each process that is needed 
to produce the functional unit, y. In the normalized system, 
the technology matrix is normalized so that each process 
gives one-unit of output. This allows the outputs to be re¬ 
moved from matrix and the inputs represented as positive 
entries (12). Each column of the normalized matrix repre¬ 
sents the inputs to produce one-unit of output of each proc¬ 
ess. Therefore the x represents the total output of each proc¬ 
ess. Any analysis should produce the same environmental 
impacts, but the outputs differ as in (13). 

The use of A and A by different authors probably results 
from different backgrounds. The matrix approach to LCA 
was originally presented in the un-normalized form (Heijungs 
1994, Heijungs and Suh 2002). IOA data is often compiled 
in an un-normalized form, but analysis is performed using 
normalized data (cf. Leontief 1970, Miller and Blair 1985, 
United Nations 1999). An advantage of using the normal¬ 
ized form is that it allows direct application of many meth¬ 
ods already established in IOA. 
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Goal, Scope and Background. To strengthen the evaluative power 
of LCA, life cycle interpretation should be further developed. A 
previous contribution (Heijungs & Kleijn 2001) elaborated five 
examples of concrete methods within the subset of numerical ap¬ 
proaches towards interpretation. These methods were: contribu¬ 
tion analysis, perturbation analysis, uncertainty analysis, com¬ 
parative analysis, and discernibility analysis. Developments in 
software have enabled the possibility to apply the five example 
methods to explore the much-used ecoinvent'96 database. 
Discussion of Methods. The numerical approaches implemented 
in this study include contribution analysis, perturbation analysis, 
uncertainty analysis, comparative analysis, discernibility analysis 
and the newly developed key issue analysis. The data used comes 
from a very large process database: ecoinvent'96, containing 1163 
processes, 1181 economic flows and 571 environmental flows. 

Conclusions. Results are twofold: they serve as a benchmark to 
the usefulness and feasibility of these numerical approaches, and 


they shed light on the question of stability and structure in an 
often-used large system of interconnected processes. Most of the 
approaches perform quite well: computation time on a moderate 
PC is between a few seconds and a few minutes. Only Monte 
Carlo analyses may require much longer, but even then it appears 
that most questions can be answered within a few hours. More¬ 
over, analytical expressions for error propagation are much faster 
than Monte Carlo analyses, while providing almost identical re¬ 
sults. Despite the fact that many processes are connected to each 
other, leading to the possibility of a very unstable system and 
very sensitive coefficients, the overall results show that most re¬ 
sults are not extremely uncertain. There are, however, some ex¬ 
ceptions to this positive message. 

Keywords: Contribution analysis; discernibility analysis; eco¬ 
invent'96; key issue analysis; life cycle interpretation; perturba¬ 
tion analysis; sensitivity analysis; uncertainty analysis 


380 


Int J LCA 12 (6) 2007 






